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Can we always get the entanglement entropy from the Kadanoff- 
Baym equations? The case of the T-matrix approximation. 
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Abstract. - We study the time-dependent transmission of entanglement entropy through an out- 
of-equilibrium model interacting device in a quantum transport set-up. The dynamics is performed 
via the Kadanoff-Baym equations within many-body perturbation theory. The double occupancy 
{hut hill) j needed to determine the entanglement entropy, is obtained from the equations of motion 
of the single-particle Green's function. A remarkable result of our calculations is that (ri.R-1-n.R4.) can 
become negative, thus not permitting to evaluate the entanglement entropy. This is a shortcoming 
of approximate, and yet conserving, many-body self-energies. Among the tested perturbation 
schemes, the T-matrix approximation stands out for two reasons: it compares well to exact results 
in the low density regime and it always provides a non-negative {hn^hm). For the second part of 
this statement, we give an analytical proof. Finally, the transmission of entanglement across the 
device is diminished by interactions but can be amplified by a current flowing through the system. 
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Originally discussed in relation to basic aspects of quan- 
tum mechanics, entanglement is nowadays seen as a key 
resource for quantum information technology [lj. Solid 
state systems are well suited to realize quantum infor- 
mation devices, for their compatibility with conventional 
electronics hardware. This has spurred a great deal of 
cross-disciplinary studies on entanglement and many-body 
physics in the condensed phase |2j. 

While ab-initio approaches to entanglement are start- 
ing to be explored [3][4], most studies so far have been for 
model systems in the ground state (see e.g. [5j|9]). How- 
ever, non-equilibrium studies are increasing in number (see 
To]- 13]. This is because knowing how entanglement 



e.g. 

is produced/transported is pivotal for controlled quantum 
information manipulations, and because it can give new 
insight into many-body systems out of equilibrium [2] . 

In this Letter we introduce a non-equilibrium Green's 
function approach to entanglement, using the time- 
dependent Kadanoff-Baym Equations (KBE) [14 15 , and 



computing a specific measure of entanglement, the entan- 
glement entropy (EE) 16] , within many-body perturba- 
tion theory (MBPT). The advantage of our approach is 
that it can be implemented ab initio |17| , if one uses re- 
liable many-body approximations (MBA:s). It is quite 
natural, at this point, to ask: When are MBA:s reliable to 




Fig. 1: (Color online) Top: model dot-leads system. Bottom: 
Time dependent local entanglement entropy £ r and double oc- 
cupancy d_R, obtained within the T-matrix and the second- 
Born approximations (TMA and BA, respectively) for a 7-site 
isolated cluster (white+red circles) with 2 particles. A per- 
turbation Vg(t) has been applied; the system parameters are 
U = 1, e = -0.5, Vo = 5, T g = 2, b g = 0.2 (see eqs.|2j4j and 
main text afterwards). 



obtain the EE ? Here, to address this issue, we choose, as 
a simple test-case model (fig. [lj, an interacting impurity 
in a quantum transport geometry. 

To anticipate the results, we see (fig. [lj, that, m some 
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time intervals, the EE may diverge, because the local dou- 
ble occupancy (used to obtain the EE) may become neg- 
ative. This non-physical behaviour is a drawback that 
MBA:s can have when computing the EE [19] . The condi- 
tions under which some standard MBA:s can or cannot be 
used to compute the EE is the main outcome of this pa- 
per. In particular, we provide an analytical proof of why, 
among all the tested MBA:s, the T-matrix approximation 
(TMA, see below) is exempt from the problem shown in 
fig. [T] 1 20 . Furthermore, our numerical results show that, 
compared to exact results in isolated clusters, the TMA 
yields the most accurate EE. As an application of our new 
method to compute the time dependent EE, in the final 
part of this Letter we will use the KBE within the TMA 
to simulate the time dependent production and transfer of 
EE in the model junction of fig. [I] and to investigate how 
the EE is affected by interactions and by a current flowing 
through such a device. Before discussing these points in 
detail, we briefly describe the model and the method. 

The model device. — We consider a single interact- 
ing impurity (at site 0) coupled to two semi-infinite, ID 
leads (figjlj • The Hamiltonian of this system is given by 



H — Hi, 



a=S.D 



[H a + V a , imp + W a (t)] + V g (t). (1) 



The impurity Hamiltonian is Hi mp = eo fiQ a + 
Uhofho^, where eo and U are the impurity on-site and 
interaction energies, respectively, and fio a = a 0a ao a , with 
a =f , I. The following three contributions describe the 
leads (source, S, and drain, D), their coupling to the im- 
purity and the effect of an external bias w a (t): 



H a = e a fk at<T - v(a\ a a a la+ i^ + H.c^j (2) 

a i a =l 

imp = -V^2a\ aa a Q . a + H.c. (3) 

(7 

OO 

w a (t) = w a {t)^2 ^2 ™ i °<>< 7 - 



V, 



(4) 



The terms w a (t) and e a are the time dependent bias and 
time independent on-site energy in the a-th lead, (both 
are taken constant in space in each lead), whereas the 
V and V are the hopping parameters in the leads and 
in between lead and impurity respectively. The time de- 
pendent gate voltage is applied on the third site in the 
S lead, V g (t) = v g (t) Y, a n 3s<a . U, e , V, w a (t), e a and 

v g {t) — Vbexp (~( *^ Tg ) 2 ) are given in units of V = 1. 
We take spin-up and -down electrons equal in number; 
this holds at all times (H has no spin-flip terms). 

We will also compare EE results from MBPT to exact 
ones in small clusters; for this, we consider a finite-size 
version of the system in fig. [T] keeping in each lead only 
the three sites closest to the impurity. The Hamiltonian 
of this seven-site cluster will be the same as in eqs. ([l] [2] 
[3II4I), but truncated to describe finite, three-site leads. 



Local entanglement entropy. — For any lattice site 
R in our system of eq. |T]), we consider the single-site EE 
and denote it by Sr. For a given many-body state | , and 
in terms of the local von Neumann entropy, £r is given by 
16 

£r = -Tr (p R log 2 p R ) , (5) 

where the reduced time-dependent density matrix, p R = 
Tr{fl'}P, is a contraction of the full density matrix p = 
|\E r )($|, and {R'} denotes the rest of the system, i.e. R £ 
{R>}. 

The local EE has been discussed for several systems 
studied with different methods, e.g. the Bethe Ansatz 
[6j[7], exact diagonalization [6j[8], the density-matrix- 
renormalization-group 12 , path-integral quantum field 



theory 10 , static [8] and time dependent [13] density func- 
tional theory, to mention a few. 

Here we add another entry to the list, by proposing 
a KBE approach to determine Sr. For a non-magnetic 
system of spin-1/2 fermions, the single-site EE is 

l°g 2 (Jy - d R 
- {l-nR+d R )\og 2 (l-nR+d R ) 



1G 



d R log 2 d R 



(6) 

In eq. (J6j), hr = (fiRf + Ur^) is the local single occu- 
pancy (directly obtained from the single-particle G and 
the KBE), and cLr = (fiR^fiR^) is the local double occu- 
pancy. When and how dR can be obtained via the KBE 
is discussed in the rest of the paper. 

Kadanoff-Baym equations and many-body ap- 
proximations. — The KBE determine the time evo- 
lution of the non-equilibrium, two-time, single-particle 

Green's function G(l, 2) = -i <V 7 ^(1)^(2) V where 

1 denote single-particle space/spin and time labels, r\0\t\ 
Here, T 7 orders the times t\, t 2 on the Keldysh con- 
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tour 7 (fig. [2jleft) 15 1), and the field operators are in the 
Heisenberg picture; the brackets () denote averaging over 
the initial state (or thermal equilibrium). Showing explic- 
itly only the time labels (matrix notation/multiplication 
in space and spin indexes is adopted), and specializing 
to time ti, we have (id tl — h (ii)) G (ti, t 2 ) = (5(12) + 
/ E (ti,t) G (t, t 2 ) dt. Here h is the single-particle Hamil- 
tonian and S, the kernel of the integral equation, is the 
self energy. E consists generally of two terms: i) Y>mba, 
which accounts for the interactions in the device region 
and is described within a given MBA (see below) and ii) 
E em f,, describing the effect of having semi-infinite contacts 
(leads). For non-interacting leads, the second contribu- 
tion ca n be treated exactly with an embedding procedure 

V'\ 2 g a , where g is the Green's 
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, giving E emb = J2 a 
function of the non-interacting, in general biased, uncon- 
nected lead a [25]. The initial state is the correlated [26] 
ground state (we work at zero temperature, 1/(3 — > 0), ob- 
tained by solving the Dyson equation G = G + G E[G]G 
self consistently, with (e — h — E em f,)Go = 1. In practice 
we solve these equations using an algorithm discussed in 
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Fig. 2: Contour (left) and approximate self-energies (right). 



In general, an exact determination of the many-body 
part £ is not possible, and one resorts to approxima- 
tions to include the effects of the interactions in the de- 
vice (hence the name £mba given to this part of £). 
Besides the TMA, the other two MBA:s we will consider 
are the second Born and the GW approximations (BA 
and GWA, respectively) (fig. fright)); they are all con- 
serving 



14 



i.e. quantities such as the total energy and 
momentum as well as the number of particles are con- 
stant during the time evolution. The self energy in the 
BA includes all terms up to second order. The GWA 
29 amounts to adding up all the bubble diagrams which 
give rise to the screened interaction, W — U + UPW, 
where P(12) = G(12)G(21). In this case the self energy 
is £(12) = G(12)W(12). In a spin-dependent treatment 
of the TMA |28)[30] one constructs the T by adding up 
all the ladder diagrams, T = $ — $[/T, where $(12) = 
G(12)G(12). The expression of the self energy then be- 
comes £(12) = / L/(13)G(43)T(34)L/(42)d34. For a de- 
tailed account of the different MBA:s see e.g. [28] , 

Double occupancy from the Kadanoff-Baym 
Equations: approximate vs exact results. While 
in principle G(12) gives access only to expectation values 
of one-body operators, its equations of motion (i.e., the 
KBE) permit to obtain also some quantities which involve 
expectation values of two-body operators, such as the to- 
tal energy. This route can be also used for the double 
occupancy da 31 ; for on-site (Ur ^ 0) interactions [32], 



drt = (nRfhRi) = - 



Ur 



£ M BA(13)G(31 + )d3 



RR 



(7) 

valid for both isolated and contacted (to leads) systems 
33 , 34 . Equation Q , together with eq. ^ and the 
KBE, is the proposed new method to compute the dynam- 
ical EE, in both model systems calculations and ab-initio 



treatments of real materials 32 



As a benchmark to the method, we considered £r for 
an isolated cluster (the seven-site system explicitly shown 
in fig. [TJ and computed eq. (J6| i) exactly and ii) using eq. 
Q and the KBE+MBA:s scheme. The numerical results 
for dR are presented in fig. [3] for the exact, BA, GWA and 
TMA treatments, in the low density regime. Two values of 
the interaction are considered, and dR is calculated for the 
last (i.e. the seventh) site of the cluster. As anticipated 
in fig. [TJ and also shown for a Hubbard dimcr in 19 1, 



Fig. 3: (Color online) Time-dependent double occupancy on 
the 7th site for the 7-site cluster with n = 1/7, T g = 2, Vo = 5, 
b g = 0.2. The curves correspond to exact (thick black solid), 
TMA (dashed red), BA (dotted green) and GWA (thin solid 
blue). In (a) U = 1 and in (b) (7 = 8. 



the TMA is the only scheme, among those considered, to 
produce correlation functions always positive. Moreover 
we notice the quite good agreement between the TMA and 
the exact dRis at low densities (higher densities, near half- 
filling, will be discussed later in the paper). These results 
clearly indicate that one should abandon the use of the 
BA and the GWA to compute the EE: in what follows, we 
will only consider the TMA. However, before analyzing 
the EE in the TMA, we describe in some detail why the 
TMA local double occupancies are always non-negative. 

Positiveness of T-matrix. — All the MBA:s used in 
this work are conserving. This, however, does not guaran- 
tee that other properties (e.g. spectral features [35]|36] or 
response functions) obtained with such MBA:s automati- 
cally satisfy further basic criteria. In fact, in fig. [TJ 3 we 
saw that density and pair correlations for GWA and BA 
may violate the positivity condition 



(^(1)^(2)^(2)^(1)) >0. 



(8) 



The pair correlation function is actually a rather sensitive 
measure of the quality of a many-body approximation. It 
has been known for a long time that the random-phase 
approximation (RPA) [37] in the ground state gives pair 
correlations which, at metallic densities, are strongly neg- 
ative at short distances [38] . Subsequently, negative pair 
correlations at short distances were also found within the 
self-consistent GWA 



:-!(> 



One may ask if the situation would improve if the 
response function was calculated self-consistently from 
changes in £ and thereby via the Bethe-Salpeter equa- 
tion: Such a response function fulfills all macroscopic con- 
servation laws, but does not necessarily yield pair correla- 
tion functions everywhere positive. For example, the RPA 
discussed above is the response function in the Kadanoff- 
Baym sense of the conserving Hartree approximation but 
it violates the positiveness condition in eq. (|8|) 38 . 



According to figs. [ljb,c) and [5] the TMA, among the 
examined MBA:s, is the only one giving always positive 
correlation functions (and, in many cases in good agree- 
ment with exact cluster results). This is not due to for- 
tuitously chosen parameter values: we will now give proof 
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that the pair correlation from the TMA is manifestly pos- 
itive, starting with the equilibrium case. 

In the ground state, for the TMA, the pair correlation 
function is a simple contraction of the T-matrix itself, 



(nR^-n R i) = -iT RR (t,t + ) = -iT RR (t,t), 



and its positiveness is a consequence of the positiveness of 
the T-matrix spectral function. Here > (<) refers to the 
electron (hole) part. The basic building block in the TMA 
is $R R >(t - f) = -iG RR ,(t - t')G RR >(t - t'). In Fourier 
space we have 



*fl«'(c) 



C RR/ (e')de' 



e — e' + irj sgn(e' — 2/i) ' 
where the spectral function C RR i is given by 



(10) 



CW(e) 



f A RR ,(e')A RR ,(e-e')de' 



Thus, C RR '(e) is positive (negative) definite for e < 2/i 
(e > 2/i). The Dyson equation for the T-matrix is T(e) = 
$(e) — $(e)C/T(e) (here $, T, U are matrices in one-particle 
indices [39], and U RR < — U R 8 RR >). This gives 



(1 + $ r U)f$ = $^(1 - UT Q ), 



and 



T(e) - ft( e ) = [1 - fr(e)U][$(e) - &{e)][l - UT a {e)], 

where we have used the identity (1 — UT)(1 + $£/) = 1, 
and a, r have their usual meaning. The T-matrix has thus 
the spectral decomposition 



T(e) 



D{e')de' 



e — e' + irj sgn(e' — 2/i) 



(11) 



where D(e) = [l-f r (e)U}C(e)[l-UT Q (e)]. Consequently, 
D(e) is positive definite for e < 2/i, and 

(n m n Rl ) = -iT RR (t,t+) = D RR {e')de' > 0. (12) 

J — oo 

The above result remains valid for a general two-body 
interaction u(ri — T2)5(ti — i2), provided we use a sym- 
metrised TMA which includes both direct and exchange 
ladder diagrams. In this case, the basic entities are 

$(rio- x , ti, r 2 , o- 2 , ti; r 3 cr 3 , t 3 , r 4 , er 4 , t 3 ) 

= (r 1 a 1 , r 2 , cr 2 |$(ti, t 3 )|r 3 CT 3 , r 4 , cr 4 ) 

= -i{G(13)G(24) - G(14)G(23)} t2=ti=M4=t3=t , , 

and T(t, t'), which can be considered as matrices in a com- 
plete set of two-particle labels and depending on two times 
t,t' |2l] . Specializing to the ground state, and in Fourier 
space, 



f(e)=$(e)-i$(e)f(°)f(e). 



(13) 




(9) Fig. 4: Correlation part of the symmetrized TMA self energy. 



Here,r(°)(12;34) = u(12) [<5(13)<5(24) - 5(14)5(23)] is a 
symmetrized interaction. Finally, for the self energy, 

(SG)(r 1; h, CTl ;r 2 , < 2 , <r 2 ) =- l / U (13)(T y n(3)V'(l)^ t (2))d3 



(14) 



u(13)T(13 + ;23 ++ )(i3. 



In eq. ( 14 ) , the first equality follows from the equation 



of motion for G, and the second defines the symmetrised 
TMA. The correlation part of S is shown diagrammati- 
cally in fig. [4] Assuming the pair interaction invertible, 
the equal-time pair correlation is given by the spin con- 
traction of 



{^(1)^(2)^(2)^(1)) = -i{U\f(t,t+)\12), 



(15) 



The very construction of the TMA in terms of multiple 
pair scatterings suggests that the individual spin compo- 
nents in eq. (15) can be interpreted as approximations 



to spin-decomposed correlations. Thus, the close corre- 
spondence between the contracted T-matrix and the pair 
correlations remains also with a general interaction, cf. 
cqs. (15 In Fourier space, the matrix $ has a spectral 
decomposition as in eq. (10), where in this case 



(12|G(e)|34) = / [A 13 {t')A 2i {e-e')-A u {e')A 23 {e-t')]de', 



and is positive de finit e for e < 2/i. From the Dyson equa- 
tion for T (cf. eq. 13 ) we then obtain 



f( £ )-ft( e ) = [l-if( e )fW][l»( e )-|.t( e )][l-if( )ft( e )]. 

Thus T has a spectral decomposition similar to that of <t, 
and its spectral function D is positive definite for e < 2/i. 
The correlation function can now be expressed in terms of 
the diagonal element of D, 



(^(1)^(2)^(2)^(1)) = -t<12|f(t,t + )|12) 



2/i 



(12|T>(e')|12)de' > 0. 



(16) 



Thus, the spin-decomposed pair correlations are positive 
and thereby also their spin averages. 

Such result remains valid in the presence of an external 
field for t > 0. The proof is similar to the proof that mani- 
festly positive =FzS^ give manifestly positive =FiG§= via the 
KBE. To simplify the analysis we stay at zero temperature 
and deform the contour in fig. [2](left) into one which runs 
from — oo to oo and back again (the equivalence follows 
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density time 



Fig. 5: (Color online) Groundstate (a) and time dependent (b) 
entanglement for the 7-site cluster on the 4th (a) and 7th (b) 
sites, with n = 1/7 for U = 1,8. In (b), T a = 2, V = 5, 
b g = 0.2 and eo = —U/2. The curves correspond to exact 
(solid) and TMA (triangles, dots). 




Fig. 6: (Color online) Time dependent entanglement on 7th 
site. V = 5, b g = 0.2, U = 1,8 and e = -U/2. In (a) the 
leads are unbiased with T g = 2 and in (b) they are biased, 
Ws/d = ±V b sm(^t), where V b = 1, T b = 20 and T g = 22. 
The green dashed line corresponds to Vb = — 1. 



from the analyticity properties of greater/lesser functions 
when Re t < and thus before any external field is ap- 
plied). With this contour, the T-matrix equation becomes 



T = $ - -|>f (°)t, 

where now the matrix labels include also time, 

(t,ri,a 1 ,r 2 ,(J2\X\t',r' 1 ,a' 1 ,r' 2 ,a' 2 ) = 
X(t 1 ,T 1 a 1 ,r 2 ,a 2 ;t' 1 ,r' l ,a' 1 ,r' 2 ,a' 2 ). 



(17) 



This leads to f< = [1 - \ f T f (°)]$<[1 - if(°)f a ]. Any 
positive-definite G < gives a positive-definite $, and be- 
cause T r is the Hermitian conjugate of T a , T < is mani- 
festly positive, and produces positive pair correlations [40] . 
As a side remark, we could add that the BA, GWA, and 
TMA all give manifestly positive spectral densities in equi- 
librium and positive ±G^ out of equilibrium, a condition 
which not all conserving schemes obey. 

Entanglement propagation across a model de- 
vice. — Since double occupancies cannot be negative in 
the TMA, we now use this approximation to investigate 
the production, offline-transmission and monitoring of EE, 
following the onset of an external gate voltage V g (t) (fig. 
[I). In the regimes considered, correlation damping and 



metastable states play an insignificant role 27 28]. 

In the low density regime, the TMA works well also 
quantitatively |27| . Thus, a suitable regime for our sim- 
ulations is when the density impurity level in the cluster 
geometry of fig. [I] both in the ground state and during the 
dynamics, has low occupations. In that regime, the differ- 
ence between EE results obtained with small and large U 
values was found to be negligible. 

This outcome can be understood looking at fig. [5^,, 
where we show, as function of the impurity density, the 
ground state EE for different U:s at the central site of our 
7-site cluster. To produce the curves in fig. [5^, the density 
is varied continuously by tuning the on-site energy. We see 
that the dependence of the EE on the interaction strength 
is only significant in the half-filling regime 41 . This re- 



mark readily applies to the time-dependent case, since we 



are computing the EE dynamics within an adiabatic local 
density approximation [13] which uses eq. (JfjJ. 

We then consider an impurity density near half-filling. 
The inherent time dependent results are shown in fig. 
[5]d, where it is immediately apparent that, as expected 
from the entanglement diagram of fig. [5^,, electron cor- 
relations considerably affect transmission. One also notes 
that while becoming quantitatively inaccurate, the TMA 
reproduces in a qualitatively correct way the exact EE. 

The effect of an electric current. — According to 
the EE diagram in fig. [5^,, if the occupation at the 4th site 
in the cluster is close to half filling, the EE is close to being 
maximal. Hence, changes in density at that site (except 
for variations limited to the small region around n = 1, 
where the EE may have a local minimum) can only reduce 
the EE. The other sites in the cluster exhibit the same 
qualitative behavior. Thus, when the isolated 7-cluster 
device is close to half filling in the ground state (as in fig. 
[5|d), a voltage pulse on the left of the device will produce a 
local reduction of the EE. At a later time, this will induce 
a change of density to the right of the device, which also 
causes a local decrease of EE. This is the behavior observed 
in fig. [5]d, i.e. the left-to-right transmission of a local 
depression in the EE, which is sensitive to the strength of 
the interaction U at the impurity site. 

For a contacted device, one has an additional handle, 
the electric current, to choose the most suitable initial 
point (interval) in the EE diagram, and thus to tune the 
"EE background" on which a voltage pulse can be super- 
imposed. This is illustrated in fig. [6] In the left panel (a), 
we display the results for an unbiased system (no current) 
at half-filling. Here, after the pulse is applied, the mech- 
anism for EE propagation is the same as in fig. [5]d, and 
we observe the transmission of a local depression in the 
EE. To provide a simple illustration of the effect a cur- 
rent may have, we applied a bias and evolved the system 
(fig. ]6]j) semi-adiabatically to approach the steady-state 
current and EE, at t ps 20, before applying a pulse V g (t). 
In this case, the transmitted EE pulse increases in ampli- 
tude from its steady state value. For U = 1, if the bias is 
reversed the pulse arrives faster, increases in size and gets 
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broadened. This effect is, though not expected, of simple 
single-particle nature. 

Conclusions. — We presented a new method to ob- 
tain the time-dependent local entanglement entropy, via 
the Kadanoff-Baym equations and within MBPT. The 
method can be implemented ab-initio, i.e. for realistic ma- 
terials. We showed that for some commonly used MBA:s 
(specifically, the BA and the GWA), the resulting double 
occupancy may become negative, i.e. leading to a diverg- 
ing entanglement entropy. We provided a proof that this 
is never the case for the TMA: the latter always gives 
positive double occupancy and, for finite systems, good 
agreement with exact numerical results in the low density 
regime. Finally, within the KBE+TMA scheme, using a 
model interacting device in a quantum transport setup, we 
illustrated the role played by the current in altering the 
modality of entanglement transmission across the device. 
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